#Figure 2: Mean Estimates of Civilian Control for Civilianized Regimes
id<-seq(1:nrow(data))
means_s<-matrix(ncol=64,nrow=nrow(theta_s_draws))
for(jj in 1:64){
x<-theta_s_draws[,id[drift_yrs==jj & drift_ind==1 & !(data$gwf_casename %in% censored_regs)]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=64,nrow=nrow(theta_d_draws))
for(jj in 1:64){
x<-theta_d_draws[,id[drift_yrs==jj & drift_ind==1 & !(data$gwf_casename %in% censored_regs)]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:64)
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
#pdf(file="age_civilianized.pdf",height=6,width=7)
plot(NULL,# create empty plot
xlim = c(1,65), # set xlim by guessing
ylim = c(-1,2), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
abline(lm(age_data$mean~age_data$age),col="#74add1")
abline(lm(age_data$mean_s~age_data$age),col="gray60")
for (ii in 0:65){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray90"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:65){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,65,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Civilianized Regime Age \n (Dropping Censored Observations)", line =4.25,cex=1) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
#dev.off()
means_s<-matrix(ncol=100,nrow=nrow(theta_s_draws))
for(jj in 1:100){
x<-theta_s_draws[,id[drift_yrs_uncapped==jj & drift_ind==1]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=100,nrow=nrow(theta_d_draws))
for(jj in 1:100){
x<-theta_d_draws[,id[drift_yrs_uncapped==jj & drift_ind==1]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:100)
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
plot(NULL,# create empty plot
xlim = c(1,100), # set xlim by guessing
ylim = c(-1,2), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
abline(lm(age_data$mean~age_data$age),col="#74add1")
abline(lm(age_data$mean_s~age_data$age),col="gray60")
for (ii in 0:100){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray80"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:100){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,100,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Civilianized Regime Age \n (Including Censored Observations)", line =4.25,cex=1) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
dev.off()
rm(list=ls())
library(dplyr)
set.seed(1989)
#Load and prep data
setwd('~/Dropbox/ctrl_isq_final/replication/appendix/irt_models')
load('data_prepped_drop_censored.RData')
load("static_output_drop_censored.RData")
theta_s_draws<-output$x
load("drift_output_drop_censored.RData")
theta_d_draws<-t(apply(output$x,1,scale))
id<-seq(1:nrow(data))
means_s<-matrix(ncol=max(drift_yrs_uncapped),nrow=nrow(theta_s_draws))
for(jj in 1:max(drift_yrs_uncapped)){
x<-theta_s_draws[,id[drift_yrs_uncapped==jj & drift_ind==1]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=max(drift_yrs_uncapped),nrow=nrow(theta_d_draws))
for(jj in 1:max(drift_yrs_uncapped)){
x<-theta_d_draws[,id[drift_yrs_uncapped==jj & drift_ind==1]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:max(drift_yrs_uncapped))
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
pdf(file="age_civilianized_drop.pdf",height=6,width=7)
par(
family = "sans",
oma = c(0,0,0,0),
mar = c(4,4,2,2),
mfrow= c(1,1)
)
plot(NULL,# create empty plot
xlim = c(1,65), # set xlim by guessing
ylim = c(-1,4), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
for (ii in 0:65){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray90"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:65){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,65,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Years of Drift", line =2.5,cex=1.25) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
dev.off()
rm(list=ls())
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(stats4)
library(MASS)
set.seed(1989)
#Load and prep data
load("/Users/kenwick/Dropbox/ctrl_isq_final/replication/irt_modeling/data_prepped.RData")
setwd('~/Dropbox/ctrl_isq_final/replication/output/')
load("static_output.RData")
output_s<-output
load("drift_output.RData")
#Extract model estiamtes and append to data object
theta_s_draws<-output_s$x
theta_d_draws<-t(apply(output$x,1,scale))
id<-seq(1:nrow(data))
means_s<-matrix(ncol=max(observed_drift_yrs),nrow=nrow(theta_s_draws))
for(jj in 1:max(observed_drift_yrs)){
x<-theta_s_draws[,id[observed_drift_yrs==jj & drift_ind==1]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=max(observed_drift_yrs),nrow=nrow(theta_d_draws))
for(jj in 1:max(observed_drift_yrs)){
x<-theta_d_draws[,id[observed_drift_yrs==jj & drift_ind==1]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:max(observed_drift_yrs))
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
pdf(file="age_civilianized_memo.pdf",height=6,width=18)
par(
family = "sans",
oma = c(0,0,0,0),
mar = c(7,7,4,4),
mfrow= c(1,3)
)
plot(NULL,# create empty plot
xlim = c(1,65), # set xlim by guessing
ylim = c(-1,2), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
abline(lm(age_data$mean~age_data$age),col="#74add1")
abline(lm(age_data$mean_s~age_data$age),col="gray60")
for (ii in 0:65){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray90"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:65){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,65,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Drift Years \n (Including Censored Observations)", line =4.25,cex=1) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
#dev.off()
censored_regs<-data$gwf_casename[data$year==1946 & data$gwf_duration>10]
#Figure 2: Mean Estimates of Civilian Control for Civilianized Regimes
id<-seq(1:nrow(data))
means_s<-matrix(ncol=64,nrow=nrow(theta_s_draws))
for(jj in 1:64){
x<-theta_s_draws[,id[drift_yrs==jj & drift_ind==1 & !(data$gwf_casename %in% censored_regs)]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=64,nrow=nrow(theta_d_draws))
for(jj in 1:64){
x<-theta_d_draws[,id[drift_yrs==jj & drift_ind==1 & !(data$gwf_casename %in% censored_regs)]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:64)
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
#pdf(file="age_civilianized.pdf",height=6,width=7)
plot(NULL,# create empty plot
xlim = c(1,65), # set xlim by guessing
ylim = c(-1,2), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
abline(lm(age_data$mean~age_data$age),col="#74add1")
abline(lm(age_data$mean_s~age_data$age),col="gray60")
for (ii in 0:65){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray90"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:65){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,65,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Civilianized Regime Age \n (Dropping Censored Observations)", line =4.25,cex=1) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
#dev.off()
rm(list=ls())
library(dplyr)
set.seed(1989)
#Load and prep data
setwd('~/Dropbox/ctrl_isq_final/replication/appendix/irt_models')
load('data_prepped_drop_censored.RData')
load("static_output_drop_censored.RData")
theta_s_draws<-output$x
load("drift_output_drop_censored.RData")
theta_d_draws<-t(apply(output$x,1,scale))
id<-seq(1:nrow(data))
means_s<-matrix(ncol=max(drift_yrs_uncapped),nrow=nrow(theta_s_draws))
for(jj in 1:max(drift_yrs_uncapped)){
x<-theta_s_draws[,id[drift_yrs_uncapped==jj & drift_ind==1]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=max(drift_yrs_uncapped),nrow=nrow(theta_d_draws))
for(jj in 1:max(drift_yrs_uncapped)){
x<-theta_d_draws[,id[drift_yrs_uncapped==jj & drift_ind==1]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:max(drift_yrs_uncapped))
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
pdf(file="age_civilianized_drop.pdf",height=6,width=7)
par(
family = "sans",
oma = c(0,0,0,0),
mar = c(4,4,2,2),
mfrow= c(1,1)
)
plot(NULL,# create empty plot
xlim = c(1,65), # set xlim by guessing
ylim = c(-1,4), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
for (ii in 0:65){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray90"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:65){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,65,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Years of Drift", line =2.5,cex=1.25) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
dev.off()
rm(list=ls())
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(stats4)
library(MASS)
set.seed(1989)
#Load and prep data
load("/Users/kenwick/Dropbox/ctrl_isq_final/replication/irt_modeling/data_prepped.RData")
setwd('~/Dropbox/ctrl_isq_final/replication/output/')
load("static_output.RData")
output_s<-output
load("drift_output.RData")
#Extract model estiamtes and append to data object
theta_s_draws<-output_s$x
theta_d_draws<-t(apply(output$x,1,scale))
id<-seq(1:nrow(data))
means_s<-matrix(ncol=max(observed_drift_yrs),nrow=nrow(theta_s_draws))
for(jj in 1:max(observed_drift_yrs)){
x<-theta_s_draws[,id[observed_drift_yrs==jj & drift_ind==1]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=max(observed_drift_yrs),nrow=nrow(theta_d_draws))
for(jj in 1:max(observed_drift_yrs)){
x<-theta_d_draws[,id[observed_drift_yrs==jj & drift_ind==1]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:max(observed_drift_yrs))
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
pdf(file="age_civilianized_memo.pdf",height=6,width=18)
par(
family = "sans",
oma = c(0,0,0,0),
mar = c(7,7,4,4),
mfrow= c(1,3)
)
plot(NULL,# create empty plot
xlim = c(1,65), # set xlim by guessing
ylim = c(-1,2), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
abline(lm(age_data$mean~age_data$age),col="#74add1")
abline(lm(age_data$mean_s~age_data$age),col="gray60")
for (ii in 0:65){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray90"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:65){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,65,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Drift Years \n (Including Censored Observations)", line =4.25,cex=1) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
#dev.off()
censored_regs<-data$gwf_casename[data$year==1946 & data$gwf_duration>10]
#Figure 2: Mean Estimates of Civilian Control for Civilianized Regimes
id<-seq(1:nrow(data))
means_s<-matrix(ncol=64,nrow=nrow(theta_s_draws))
for(jj in 1:64){
x<-theta_s_draws[,id[drift_yrs==jj & drift_ind==1 & !(data$gwf_casename %in% censored_regs)]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=64,nrow=nrow(theta_d_draws))
for(jj in 1:64){
x<-theta_d_draws[,id[drift_yrs==jj & drift_ind==1 & !(data$gwf_casename %in% censored_regs)]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:64)
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
#pdf(file="age_civilianized.pdf",height=6,width=7)
plot(NULL,# create empty plot
xlim = c(1,65), # set xlim by guessing
ylim = c(-1,2), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
abline(lm(age_data$mean~age_data$age),col="#74add1")
abline(lm(age_data$mean_s~age_data$age),col="gray60")
for (ii in 0:65){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray90"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:65){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,65,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Civilianized Regime Age \n (Dropping Censored Observations)", line =4.25,cex=1) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
#dev.off()
rm(list=ls())
library(dplyr)
set.seed(1989)
#Load and prep data
setwd('~/Dropbox/ctrl_isq_final/replication/appendix/irt_models')
load('data_prepped_drop_censored.RData')
load("static_output_drop_censored.RData")
theta_s_draws<-output$x
load("drift_output_drop_censored.RData")
theta_d_draws<-t(apply(output$x,1,scale))
id<-seq(1:nrow(data))
means_s<-matrix(ncol=max(drift_yrs_uncapped),nrow=nrow(theta_s_draws))
for(jj in 1:max(drift_yrs_uncapped)){
x<-theta_s_draws[,id[drift_yrs_uncapped==jj & drift_ind==1]]
means_s[,jj]<-apply(x,1,mean)
}
means<-matrix(ncol=max(drift_yrs_uncapped),nrow=nrow(theta_d_draws))
for(jj in 1:max(drift_yrs_uncapped)){
x<-theta_d_draws[,id[drift_yrs_uncapped==jj & drift_ind==1]]
means[,jj]<-apply(x,1,mean)
}
age<-seq(1:max(drift_yrs_uncapped))
age_data<-as.data.frame(age)
age_data$mean_s <-  apply(means_s, 2, mean)
age_data$mean_up_s <- apply(means_s, 2, quantile, 0.975 )
age_data$mean_lw_s <- apply(means_s, 2, quantile, 0.025 )
age_data$mean <-  apply(means, 2, mean)
age_data$mean_up <- apply(means, 2, quantile, 0.975 )
age_data$mean_lw <- apply(means, 2, quantile, 0.025 )
age_data$age_s = age_data$age-.25
age_data$age_d = age_data$age+.25
pdf(file="age_civilianized_drop.pdf",height=6,width=7)
par(
family = "sans",
oma = c(0,0,0,0),
mar = c(4,4,2,2),
mfrow= c(1,1)
)
plot(NULL,# create empty plot
xlim = c(1,65), # set xlim by guessing
ylim = c(-1,4), # set ylim by the number of variables
axes = F, xlab = NA, ylab = NA)  # turn off axes and labels
grid()
for (ii in 0:65){
lines(c(age_data$age[ii]-.25,age_data$age[ii]-.25),c(age_data$mean_up_s[ii], age_data$mean_lw_s[ii]),col="gray60", lwd=1.75)
points(age_data$age[ii]-.25,age_data$mean_s[ii],bg=c("gray90"), col="gray60",  cex=.65, pch=21, lwd=1.5)
}
for (ii in 0:65){
lines(c(age_data$age[ii]+.25,age_data$age[ii]+.25),c(age_data$mean_up[ii], age_data$mean_lw[ii]),col="#313695", lwd=1.75)
points(age_data$age[ii]+.25,age_data$mean[ii],bg=c("#74add1"), col="#313695",  cex=.65, pch=21, lwd=1.5)
}
lab=c(0,(seq(5,65,by=5)))
axis(1, at=lab, labels=lab,cex.axis=1)
axis(4,at=NULL, labels=FALSE,tick=FALSE,col="white",cex.axis=1)
axis(side=2,  las=1,cex.axis=1)
mtext(side = 1, "Years of Drift", line =2.5,cex=1.25) # label bottom axis
mtext(side = 2, "Mean Civilian Control", line=2.25, cex=1.25)# add title
text(0,2.75,"Static",pos=4,col="gray60",cex=1.5)
text(0,2.5,"Dynamic Drift",pos=4,col="#313695", cex=1.5)
box()
dev.off()
